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We describe here in detail the recently introduced methodology for simulation of structural tran- 
sitions in crystals. The applications of the new scheme are illustrated on various kinds of crystals 
and the advantages with respect to previous schemes are emphasized. The relevance of the new 
method for the problem of crystal structure prediction is also discussed. 



I. INTRODUCTION 



The constant-pressure molecular dynamics (MD) 
method of Parrinello and Rahman[l| enabled for the first 
time the study of structural phase transitions in bulk 
solid crystals by computer simulation. Starting from a 
known initial structure it allowed the identification of 
possible candidates for the new structure without any 
previous knowledge. It thus achieved predictive power, 
in particular in combination with ab-initio methods 2], 
and has been successfully applied many times to a vari- 
ety of crystalline systems (for few selected applications, 

• see Refs.]! S H S B 0). In the practical use of the 
method, however, several problems arise related mainly 
to the fact that structural transitions are often first order. 

\ This is by necessity the case when the symmetries of the 
two crystal structures are not in a group-subgroup rela- 
tion. Experimentally, such transitions proceed via nucle- 
ation of the new phase, which often starts on the surface 
or on structural defects. For simulations of crystals, how- 
ever, periodic boundary conditions that eliminate surface 
are commonly used. The systems used in simulations 

\ are typically relatively small and therefore contain no 
structural defects. The simulation setup therefore sup- 
presses the possibilities for a heterogeneous nucleation of 
the new phase. As a consequence, the transformation 
of the system proceeds in a collective way, involving all 
atoms, which results in a high barrier. This might cause 
a metastability of the initial phase far beyond the ther- 

\ modynamic transition point and large hysteretic effects 
are frequently observed. Should there be a substantial 
difference in volume between the phases, increasing pres- 
sure favors the new phase with a smaller volume, due to 
the contribution of the PV term in the Gibbs potential. 
In order to observe the transition within the accessible 
simulation time one often has to overpressurize the sys- 
tem close to the point of mechanical instability Un- 
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der such conditions one or more intermediate phases may 
be skipped j2L I lln]. which reduces the predictive power of 
the method. In other cases, where the volume difference 
between the phases is not so pronounced, even overpres- 
surization might not help to force the transition to occur. 
For the above reasons, there is still a need for developing 
specific methods aimed at simulating structural transi- 
tions in crystals, given the great theoretical and prac- 
tical relevance of the closely related problem of crystal 
structure prediction. 

In this paper we review recent progress in this field re- 
lated to the application of the new approach of Laio and 
Parrinello, called metadynamics [l^ . Rather than giving 
a ready-to- use recipe we would like to highlight the possi- 
bilities offered by the new methodology which allows the 
design of suitable algorithms for different kinds of sys- 
tems. In section^ we review the generic algorithm de- 
veloped in Rcf. 13] in which metadynamics is performed 
using the simulation box as order parameter. The use of 
the method will be illustrated with examples of zeolite 
and benzene crystals. In section lTlI Al we discuss a differ- 
ent variant of the method (Ref. 14]) suitable for systems 
where an internal order parameter has to be used instead 
of the simulation box. This case is illustrated by the ex- 
ample of graphite-to-diamond conversion. Finally, in the 
last part we draw some conclusions and suggest possible 
directions for further development. 

II. METADYNAMICS USING A SIMULATION 
BOX AS ORDER PARAMETER 

Since the Parrinello-Rahman method represents a 
generic constant-pressure MD simulation method, the 
problems related to its application to structural transi- 
tions originate in the lack of efficiency of standard MD 
in crossing high barriers. During a standard equilibrium 
simulation the system explores only a small part of its 
free energy surface, corresponding to thermal fluctua- 
tions around a locally stable minimum. Consequently, 
a spontaneous passage to another minimum separated 
by a barrier substantially larger than the thermal energy 
/csT is extremely unlikely on the time scale of a typical 
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MD simulation, which is of the order of ns in the case of 
classical force fields and ps in the case of ab-initio meth- 
ods. 

A new general approach for escaping free-energy min- 
ima and systematic exploration of free-energy surfaces 
has recently been developed by Laio and Parrinello[l^. 
It has been called metadynamics and is capable of dra- 
matically speeding up the simulation of activated pro- 
cesses involving barrier crossing [Tsl ITA IT^ . including 
first-order phase transitions. The general algorithm has 
been adapted by Martohak, Laio and Parrinello for 
simulating structural transitions in crystals. We note 
that this method represents a conceptual extension of 
the idea of constant-pressure simulation. Instead of the 
latter, we perform a search for new minima by exploring 
the surface of the Gibbs potential. In order to reduce the 
complexity of the problem, an order parameter or col- 
lective variable is needed that allows us to discriminate 
between different crystalline structures. A natural one 
for a crystal is the structure factor ^(k) which provides 
a unique fingerprint of the spatial arrangement of the 
atoms in a periodic lattice. The structure factor, how- 
ever, is not a convenient order parameter to use for our 
purpose, because it has high dimensionality (in principle 
infinite). Instead, we follow in this point the approach of 
Parrinello and Rahman and use as order parameter di- 
rectly the three supercell edges a, b, c, arranged as a 3 x 3 
matrix h = (a, b,c). For relatively small systems, where 
the creation of defects is too costly, the box matrix h is 
likely to be simply related to the unit cell u via relation 
h = um, where m is an integer matrix. The matrix h 
can therefore distinguish between different unit cells and 
crystal structures. Out of the 9 independent degrees of 
freedom of the matrix h only 6 determine the shape of 
the box, while the remaining 3 are related to the global 
rotation of the box. The latter is irrelevant and only com- 
plicates the analysis of the results. Therefore it is conve- 
nient to freeze it, thus reducing the number of degrees of 
freedom to 6. In Ref.^^ we followed the idea of Nose and 
Klein Il9| and used a symmetric matrix h. An equally 
good choice is to constrain h to an upper triangular form 
and here we follow this way. Since our aim is to simulate 
a phase transition at a pressure P and a temperature T 
we will explore the Gibbs potential ^(h) = !F{h)+PV as 
a function of h where .F(h) is the Helmholtz free energy 
of the system at fixed box and V = det(h) is the volume 
of the box. The six non-zero components of h can be 
conveniently arranged as h = (/in, /122, hi2, his, ^23) 
and represent a 6-dimensional order parameter. This dis- 
tinguishes the different local minima of Q corresponding 
at pressure P to stable or metastable crystal structures. 

The metadynamics requires calculation of derivatives 
of the free energy with respect to the order parameter. 
In our case such a derivative has a simple form 

where p is the internal pressure tensor, which can be 



easily evaluated in MD or Monte Carlo simulations at 
constant h from the averaged microscopic virial tensor 

The order parameter is now evolved according to a 
steepest-descent-like discrete evolution with a stepping 
parameter Sh (metadynamics) 

h-^=h* + J/.|^. (2) 

Here, the driving force = derived from a 

history-dependent Gibbs potential 0* where a Gaussian 
has been added to G (h) at every point h* already visited 
in order to discourage it from being visited again. Hence 
we have 

g*(h) =g(h) + ^VFe-^'5?s^ (3) 
t'<t 

and the force is thus a sum of a thcrmodynamical 
driving force F — and the term Fg coming from a 
potential constructed as a superposition of Gaussians. As 
time proceeds the history-dependent term in Eq.Q fills 
the initial well of the free-energy surface and the system 
is driven out of the local minimum. 

The metadynamics algorithm described above can be 
implemented as follows (Fig^J. We start from an equi- 
librated box h containing the initial structure at a given 
pressure P and temperature T and evaluate the pressure 
tensor p in a constant h MD run long enough to allow re- 
laxation to equilibrium and sufficient averaging of p. The 
box h is then updated using the forces and metady- 
namics equations (j2l3(l to a new value h . After the box 
is modified the particle positions are rescaled in order 
to fit into the new box using the relation r = h h~^r. 
In the case of molecular crystals we only scale the cen- 
ters of mass of the molecules but not the intramolecu- 
lar degrees of freedom. As the initial free energy well is 
gradually filled the box undergoes a sequence of progres- 
sively larger deformations until a transition occurs and 
the system enters into the basin of attraction of a new 
state. This can be detected by monitoring the structure 
factor S'(k) and is often apparent also on a visual inspec- 
tion of the atomic configuration. At this stage one can 
switch off the Gaussian term, so that the metadynamics 
becomes purely steepest-descent-like and drives the sys- 
tem towards the equilibrium state for the new structure. 
In this equilibrium state the pressure will be equal to P. 
Once the new structure is characterized one can switch 
the Gaussians on again, thus filling the new minimum, 
and move to other minima, if available. 

The important point for the success of the method is 
the judicious choice of the parameters W and 6h. In 
principle, these depend on the t/(h) landscape. The pa- 
rameter Sh determines the resolution in h and should be 
smaller than the typical size of the well. A simple way 
to estimate the order of magnitude of Sh is to perform 
a short Parrinello-Rahman simulation and compute the 
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start from an initial structure 



equilibrate the system 



average the pressure 
tensor p 




switch off the gaussians 
(set W = 0) 



equilibrate the system 



average the pressure 
tensor p 



calculate a new box h' 
according to Eq. (2) 



calculate a new box h' 
according to Eq. (2) 



rescale the positions 
of the particles to h' 



rescale the positions 
of the particles to h' 




stop 



FIG. 1: Flowchart diagram of the metadynamics algorithm 
for seeking new crystal structures. 

fluctuations of the box h. However, we note that a very 
small value of Sh is not likely to be useful. Since the 
volume of the 6-dimensional Gaussian is proportional to 
the sixth power of its size, a small Sh means that many 
Gaussians would be needed to fill the well, resulting in a 
long run. In order to achieve the necessary energy reso- 
lution, W should be chosen as a fraction of the relevant 
energy barriers. These are, however, not known at the 
beginning, since one does not know the exact mechanism 
of the transition to the new phase. A practical guideline 
for the choice of W and Sh can be based on the require- 
ment that the Gaussians should not be too "sharp", or 
in other words, for an optimal filling the curvature of the 
Gaussians should be smaller than that of the well (see 
Ref. jl^j for a more general discussion of the choice of W 
and Sh). This leads to the condition < K where K is 
the smallest eigenvalue of the Q(h) Hessian at the min- 
imum ho. For a cubic system we can estimate K from 
the approximate expansion of t?(h) around ho 



g(Ah)«g(ho) + iFc(^)2 



(4) 



of the elastic constants. This provides an estimate K 
Lc and a condition 

W 



Sh^ 



< Lc 



(5) 



where L is the cell edge and c is of the order of magnitude 



we stress here that the right-hand side of the inequality 
contains the box size L and therefore the choice of the 
parameters W and Sh is dependent on the size of the 
simulation box. In the examples that follow the criterion 
(|5|) was found to provide a good working procedure. 

We now discuss several aspects of the algorithm. First 
we note an iinportant difference with respect to the gen- 
eral method |l2l |. In principle, metadynamics may al- 
low the free energy profile to be recovered. As shown in 
Ref. [l^ ■ the sum of the Gaussians in Eq.© converges 
to — t/(h) up to an additive constant, provided W and 
Sh are properly chosen, the system is confined to a finite 
region of the order parameter space and several crossings 
between the minima occur. Such information would, of 
course, be very valuable, since it would allow an accu- 
rate determination of the transition pressure. The lat- 
ter is, in general, difficult to calculate once the anhar- 
monic effects start to play a role and one cannot ap- 
ply the quasiharmonic or self-consistent quasiharmonic 
approximation. The recovering of the free energy, how- 
ever, implicitly assumes (similarly to thermodynamic in- 
tegration) that the system always evolves in a reversible 
way, and the free energy never decreases discontinuously. 
This is generally not guaranteed in the case of solid-solid 
transitions when box h is used as an order parameter. 
In this case the metadynamics should be effectively re- 
garded as an "engine" which by inducing a progressively 
larger deformation of the initial structure guarantees at 
some point the transition to a new one. Whether such a 
transition proceeds in a reversible or irreversible way, as 
well as the particular mechanism involved, is clearly very 
much system dependent. It can proceed via the creation 
of defects, such as grain boundaries or stacking faults, 
or in a collective way, as when a mechanical instability 
is reached by increasing the pressure in the Parrinello- 
Rahman method. There is, however, a substantial differ- 
ence related to the symmetry. In the constant-pressure 
Parrinello-Rahman simulation the internal stress tensor 
is on average equal to a prescribed value, which is usu- 
ally chosen as the hydrostatic pressure P. Such compres- 
sion preserves the initial symmetry of the system and un- 
less the system is close to an instability point, the sym- 
metry is only slightly broken by instantaneous thermal 
fluctuations. Metadynamics, however, is not a constant- 
pressure simulation but rather a method aimed at ex- 
ploring the free-energy surface. Starting from an initial 
minimum, the exploration consists of applying various 
symmetry-breaking deformations of the initial structure, 
which clearly facilitates reaching a structure with a dif- 
ferent symmetry. The internal stress tensor during the 
exploration may become substantially anisotropic and its 
fluctuations around the prescribed (hydrostatic) value P 
are much stronger than the thermal ones. For this rea- 
son, the role of P does not appear to be critical. As 
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long as a given structure is at least metastable at pres- 
sure P, it can in principle be found by our algorithm; 
as discussed at the beginning of the previous section, the 
phases are usually metastable in a broad pressure interval 
around the equilibrium transition pressure Pc- However, 
the value of P affects the number of steps necessary to 
reach the transition. If its value is too low compared to 
Pc, a lot of Gaussians may be needed to fill the initial 
well before a transition is observed. 

A second remark should be made concerning the pos- 
sibility that the system might change box h without re- 
ally changing the structure. This usually proceeds via 
plane sliding and is a manifestation of the fact that one 
structure can be described by many equivalent choices 
of box h, known as modular invariancepl}. This phe- 
nomenon cannot be avoided completely as it is related to 
the very nature of the order parameter h. Its frequent 
occurrence, however, together with a failure of the algo- 
rithm to find a new structure, might be a manifestation 
of either too wrong a pressure or the fact that box h is 
not a good order parameter for a given system (see sec- 
tion lIII|l . According to our experience the choice of a not 
too small value of 5h, inducing substantial volume fluc- 
tuations, helps to suppress the transitions to the same 
structure, as these obviously conserve the volume. Too 
large a value might, however, cause a transition to an 
amorphous structure; a suitable compromise has to be 
found by trial and error. 

Our original application to an atomic system, silicon 
crystal, is described in the paper 13|; the method worked 
very well in that case. To assess its applicability to a 
wider class of systems we also performed a case study 
of two other crystals of a rather different kind. These 
are presented in the following subsections. Zeolite is an 
example of a crystal having an extended network while 
benzene represents an organic molecular crystal. 



A. Application to zeolite 

As our first example we provide here a summary of 
the main results obtained by the application of the above 
method to reconstructive transitions on a complex frame- 
work structure like a zeolite|22|. 

Zeolites are tectosilicates, with a structure formed by 
corner sharing Si04 or AIO4 tetrahedra and characterized 
by having cages or channels able to accommodate alkaline 
or earth-alkaline cations and small molecules (generally 
water) "2?, '24*1 . We have focused on the Li-ABW ze- 
olite (Li[AlSi04]-H20), first synthesized by Barrer and 
White in 1951 [23 ■ The framework structure has an 
orthorhombic symmetry [26l | and is formed by directly 
connected hexagonal rings sheets in the be plane. Wa- 
ter molecules and Li cations are located inside monodi- 
mensional 8 membered channels developing along the c 
axis. If the temperature is raised in dry environment, the 
Li-ABW zeolite undergoes two phase transitions. A first 
displacive transition occurs to the framework when dehy- 



dration starts, leading to a structure known as anhydrous 
Li-ABW. At a temperature around 650° C a second recon- 
structive transition drives the system to a new structure, 
the 7-eucryptite. This transition is accompanied by a 
symmetry change, leading to a monoclinic crystal. The 
framework, formed by 8- and 4-membered rings of tetra- 
hedral units in the anhydrous Li-ABW, transforms in to 
a 6-membered rings structure. Because of the high sta- 
bility of the Al-0 and Si-0 bonds a Parrinello-Rahman 
approach is not able to reproduce the reconstructive tran- 
sition; in order to provide enough thermal energy to ob- 
serve a spontaneous breaking of some bonds, the temper- 
ature has to be increased so much that only a collapse of 
the structure is observed. 

The simulation of the anhydrous Li-ABW — > 7- 
eucryptite transition has been carried out within the 
metadynamics approach employing a classical potential, 
using the form and the parameters proposed in literature 
by Zirl and Garofalini Details about the structure 
of the potential and its adaptation for our purposes can 
be found in Ref . p3| . 

The metadynamics run was performed at external 
pressure P = for a system consisting of 56 atoms, us- 
ing the parameters for the Gaussian term equal to W 
— 9.94kcal/mol and 5h — 0.4 A. The total simulation 
time for each NVT run was 10 ps with an equilibration 
of 2.5 ps. During the first 67 steps of metadynamics 
the deformation of the simulation cell does not lead to 
bond breaking, but only to displacive rearrangement of 
the framework. After step 67 a reconstructive event is 
recognized by a sudden drop of the strongest diffraction 
peaks of the original structure (see Fig[21). The transition 
from anhydrous Li-ABW to 7-eucryptite is completed in 
three metadynamics steps and it involves only processes 
of bond switching between some Al-0 bonds (Fig|21l. In 
particular, half of the Al atoms of the simulation cell 
move in a concerted way, leading to four bond switch, 
such that the Li-ABW tetrahedral units are broken. Af- 
ter few picoseconds, a second bond switch is observed 
for each of the Al atoms already involved in the first 
step, leading to the formation of new tetrahedral units 
with the topology of the eucryptite. While the Al atoms 
involved in the bond breaking move significantly, explor- 
ing three main positions, the other atoms at the center 
of the unbroken tetrahedral units move very little from 
their starting crystallographic positions, in spite of the 
fact that each tetrahedron undergoes a rigid rotation in 
order to complete the transition. The Li atoms also par- 
ticipate in the transformation with a displacement from 
their starting position near to the Al involved in the bond 
switch. The transition pathway observed with metady- 
namics is in good agreement with the experimental obser- 
vations reported in Ref. . A more detailed description 
of the calculations and the transition path are given in 
Ref. [22. 

The above results demonstrate the ability of metady- 
namics to uncover the detailed microscopic mechanism 
of a phase transition in a complex crystal structure. An- 
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FIG. 2: Zeolite: time evolution of cell lengths, angles and 
selected strong peaks of the structure factor during the meta- 
dynamics simulation. 



other attractive feature of the procedure we outlined in 
the previous section is that it allows a systematic ex- 
ploration of the different polymorphs in which a given 
system can exist. Starting from an initial structure, sev- 
eral metadynamics runs can be performed at different 
conditions of pressure and temperature until a transition 
is observed. Any classical potential providing a quali- 
tative description for the system can be employed, and 
performing the metadynamics with an accurate model 
(e.g. DFT) might be a waste of computer time. In fact, 
the Gibbs free energy difference between the candidate 
structures found by metadynamics can be computed with 
the high level method once the structures are known at 
a relatively small computational price. 

For the system under investigation, we performed sev- 
eral metadynamics runs starting from the Li-ABW struc- 
ture in the pressure range between 10 and 100 kbar, ob- 
taining seven new structures characterized by a differ- 
ent connectivity of the atoms. To illustrate the excel- 
lent capability of the technique to find new polymorphs 
we report here a brief description of two of the struc- 
tures we find (these results will be described in detail in 
a separate publication H^). The phase pictured in Fig^ 
(a) is formed by ordered sheets, balanced by the pres- 
ence of Li cations in the intra-layers space. Each layer is 
formed by alternating Si-centered and Al-centered tetra- 
hedra that form four-membered rings. The structure 
reported in Fig^ (b) has the same connectivity of 7- 
eucryptite, but in this case there is no perfect alterna- 
tion of Si-centered and Al-centered tetrahedra, whereas 
there are some Al-O-Al and Si-O-Si bridges. The Al-0- 
Al bridge is normally not found in zeolitic frameworks 
(Lowenstein's rulc^29j) because it shows larger stability 
for angles around 180° which are generally not present in 
zeolites, but in conditions of high external pressure such 
configuration can be stabilized. For these two phases and 
the other five new structures found by the metadynamics 
procedure we determined the ab-initio equation of state 





FIG. 3: Zeolite: three main steps involved in the transition 
path. The blue tetrahedra are the Al-centered ones, while the 
grey tetrahedra are those that are Si-centered. The tetrahe- 
dra directly involved in the bond breaking are pictured in ball- 
and-stick representation, where the red spheres represent the 
oxygen atoms. The Li cations inside the cavities are omitted 
for clarity, (a) The anhydrous Li-ABW structure (8- and 4- 
membered rings), (b) an intermediate structure, formed dur- 
ing the first bond switching, (c) the final 7-eucryptite struc- 
ture (6-membered rings). 



and the phonon spectrum in order to estimate the en- 
tropic contribution to the free energy, thus obtaining a 
phase diagram for the system at the DFT level and with- 
out any a priori knowledge of the possible phases of the 
system 28J . 
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FIG. 4: Zeolite: two phases obtained by performing metady- 
namics at high external pressure. More details will be pro- 
vided in a forthcoming publication|28j|. 



B. Application to benzene 

The possibility of predicting all of the stable crys- 
talline structures and thus the physical and chemical 
properties for a given compound is extremely impor- 
tant for pharmaceuticals |^ and benzene is a good 
case study to prove the applicability of our technique 
to organic molecules. Benzene, though being a sim- 
ple molecule has a quite complex phase diagram. Since 
the late 60's many authors studied solid benzene and in 
the literature two different notations for its crystalline 
phases can be found. In this paper we follow the no- 
tation used by Thiery and Leger in Ref. (32]. At room 
temperature and ambient pressure benzene crystallizes 
in an orthorhombic structure, Pbca, which is stable up 
to 1.4 GPa j32|. Beyond this pressure a sluggish transi- 
tion to benzene II is observed [s^ [H, |3 ■ Determining 
the crystalline structure of this phase was a big challenge 
and only 10 years after its first observation theoretical 
calculations could determine that it is orthorhombic and 
belongs to the P43 2i2 space group Benzene II is 

stable up to 4 GPa and then it transforms into benzene 
III, which is monoclinic and belongs to the P2\/c space 
group [3^ lU 13 . Upon increasing pressure two more 
phases are observed: benzene IIP, which is stable be- 
tween 11 and 24 GPa, and benzene IV, which is stable 
at even higher pressure [s^, 13- These latter two struc- 
tures are not well characterized and there is still some 
debate on whether they are real phases or not. In partic- 
ular, benzene III' is supposed to be only a modification of 
benzene Illjs^ and benzene IV to be polymorphic [s^. 
Beyond these five, another phase has been hypothesized: 
benzene P, which should be stable at low temperatures 
(^100 K) 32] and normal pressure. 

In this work we focus our attention only on benzene 
I, II and III. These arc indeed the only three phases for 
which the structure is well understood both by experi- 
ment and by theory. In a different publication we tackle 
the problem of determining all stable structures of ben- 
zene [3^ and show that the method can indeed identify 
the structure of all experimentally proposed phases. 



All our MD simulations were done using the GRO- 
MOS96 force field 37]. The molecule is fully flexible and 
the long range electrostatic interaction is calculated with 
the Particle Mesh Ewald (PME) summation. We note 
that the intermolecular potential can reproduce the sta- 
bility of all the three phases, though not their correct 
transition pressures. At K, benzene II is less stable 
than benzene III and, at variance with the experiments, 
the common tangent construction suggested a transition 
pressure of about 1.5 GPa for both the I to II and I to 
III transformationist, ■ 

We apply metadynamics to study the transitions 
among those phases. During the metadynamics both the 
equilibration and the averaging of the internal stress ten- 
sor are performed by NVT runs of 1 ps, the external 
pressure P is always fixed at 2 GPa and the temperature 
at 300 K. For what concerns the parameters entering the 
time dependent potential, we tried a wide range of val- 
ues, however, good results are obtained only when the 
Gaussian width, is chosen between 2 and 3 A and 
the Gaussian height, between 120 and 600 kcal/mol. 
Given the elastic constants of benzene of about 6 GPa 
and the box vectors varying between about 10 and 50 A 
we notice that the working parameters fulfill the empiri- 
cal guideline given in section^ (Eq. [SJ. We experienced 
that the correct choice of 8h is extremely important. In 
fact, if i5/i < 2 A during metadynamics run we only ob- 
serve a useless plane sliding, which is indeed a much less 
expensive deformation than a genuine solid-solid phase 
transition. For molecular solid this is even more easy 
than in covalent crystals because of the weak nature of 
the intermolecular interactions. Moreover, if i5/i > 3 A 
the deformation of the cell was too large to be accom- 
modated by the molecules and amorphization is always 
observed. On the other hand, choosing a too small value 
for W causes a slow escaping from the energy basin re- 
sulting in a waste of time. 

The structure factor provides a unique fingerprint of 
each crystalline structure and can be compared with the 
experimental x-ray powder diffraction pattern. Charac- 
teristic reflection peaks of a given structure often vanish 
for other crystalline structures and by monitoring the 
intensity of the strongest reflections during the metady- 
namics loop we can easily detect when a phase transition 
occurs. For a typical metadynamics run, where the phase 
I transforms into phase III, we show in Figure the be- 
havior of the box lengths and angles (top panels) as well 
as the intensity of a selected peak of the structure factor, 
relative to the strongest reflection. 

We tried tens of starting configurations with either of 
the three phases and with different dimensions of the 
simulation cell, containing from 4 to 256 molecules. We 
observed that by increasing the cell dimension the num- 
ber of metasteps necessary to escape from the energy 
basin decreases. At the same time the formation of de- 
fects becomes easier and in the cells larger than about 
hundred molecules we sometimes observe stacking faults. 
We stress, however, that we never observe "point" de- 
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FIG. 5: Benzene; behavior of the box lengths and angles in 
a typical metadynamics run with 192 molecules in the sim- 
ulation cell (top and middle panels). In the bottom panel 
we report the intensity of the (434) reflection relative to the 
strongest peak observed in that metastep. The sudden drop 
of the intensity clearly indicates the occurrence of a phase 
transition at the metastep 22. 



fects, such as vacancies or misoriented molecules. In par- 
ticular, we could never observe a clean transition between 
benzene I and II if the simulation cell was larger than 32 
molecules. With larger cells benzene I always transforms 
either to benzene III, which is lower in energy than ben- 
zene II, or to a defective structure. In Figure we show 
the plane stacking along the c axis for benzene II and III 
and also for a typical defective configuration obtained by 
metadynamics . 

For comparison, we performed also a limited amount 
of Parrinello-Rahman simulations, starting from phase 
I. We observed in this case only a transition to phase 
III, but not to phase II. Moreover, no phase transition is 
observed at pressures lower than 10 GPa. The results of 
metadynamics simulations clearly demonstrate that the 
new method works very well also for molecular crystals. 



III. METADYNAMICS USING AN INTERNAL 
ORDER PARAMETER 

The new method described in section ^1 still suf- 
fers from some limitations in common with the original 
Parrinello-Rahman method. For instance it is less effec- 
tive for the study of phase transitions for which the pri- 
mary order parameter is an internal coordinate instead of 
the cell edges. This is the case for phase transformations 
under pressure described in terms of solid state chemical 
reactions such as the 2D ,38j and 3D 39] polymerizations 
of Ceo or the topochemical solid-state polymerizations of 
alkenes, alkynes and aromatic hydrocarbons H,E3- ^'^^ 
instance, in the 2D polymerization of Cgo the activation 
barrier for the [2-1-2] cycloaddition reaction is overcome 
by a suitable deformation of the fuUerenic cage which is 
not induced by simply decreasing the intcrcage distances 
down to the density of the 2D polymer |_41]. In the per- 
spective to address the study of phase transformations in 
this class of materials, we have extended the metadynam- 
ics scheme recently devised by lannuzzi, Laio and Par- 
rinello (ILP) ^2] to constant-pressure MD simulations 
01 . The ILP scheme can be dubbed reactive molecular 
dynamics since suitably defined reaction coordinates are 
introduced as dynamical variables. 

By combining the ideas behind the Parrinello-Rahman 
and ILP methods we have introduced a constant-pressure 
reactive MD described by a Lagrangian of the form 

1 ^ 1 
C = -J2m^{slh'hs,)-E{{s,},h) + -W,Trh'h 

1=1 



V{t,{Vc.}), 



(6) 



where the first line is the Parrinello-Rahman Lagrangian 
and the second line is the ILP Lagrangian Here, 
Si are scaled ionic coordinates, fl is the cell volume, p 
the external pressure and r]a are collective variables as in 
the ILP scheme ^3] with a fictitious kinetic energy and 
mass {Ma). A harmonic potential restrains the values of 
the collective coordinates r]a{{si},h.) close to the corre- 
sponding dynamical collective variables rja- The values 
of Ma and ka control the time scale of the evolution of 
the collective variables and are chosen according to the 
prescription given in Ref. [42| . The collective coordi- 
nates r]a{{si},h) are functions of the scaled ionic coordi- 
nates and of the cell edges and should be able to discrim- 
inate between the initial and final phases. E{{si},h) is 
the total internal energy while V{{r]a},t) is the history- 
dependent potential acting on the collective variables and 
given by 



V{t,{rja}) = J2wl[< 



(7) 



t'<t 



where W and CTq are suitably chosen parameters as de- 
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FIG. 6: Benzene: comparison of the molecular arrangement of phases 11 and III with that observed in a typical defective 
structure obtained during metadynamics. Note that benzene II has a typical abed plane stacking while benzene III has an a 
b stacking. The defective structure can be viewed either as a crystal of benzene II with two stacking faults (along the planes 
marked by a star) or as benzene III with two stacking faults (in the planes marked by a prime). 



scribed in Ref . ■ The equations of motion correspond- 
ing to the Lagrangian © are reported in Ref. Q . 



A. Application to graphite 

To demonstrate the validity of the new scheme de- 
scribed above we have simulated the conversion of car- 
bon from graphite to diamond under pressure p^ . This 
transformation can be seen as driven by an internal order 
parameter such as the corrugation of the graphitic planes 
which leads to the change of hybridization of carbon from 
sp^ to sp^. As a measure of the hybridization type of the 
carbon atoms, we have defined as collective variable the 
coordination number of the atoms of a single graphitic 
plane in the simulation cell with respect to the atoms of 
the two neighboring planes, i.e. 



with 



E E ' 

i^plane j^plane 



(8) 



(9) 



where is the distance between atoms i and j and 
d — 2.2 A. The graphite to diamond conversion has been 
already reproduced in the ab-initio Parrinello-Rahman 
molecular dynamics simulations of Ref. , although at 
a pressure of 90 GPa, six times larger than the experi- 
mental estimate (15 GPa, 9]), due to the aforementioned 
limitations of the Parrinello-Rahman method. In the 
present work graphite is described by the tight-binding 
(TB) potential of Ref. ;43] supplemented by an empiri- 
cal two-body van der Waals (vdW) interaction, neces sary 
to describe the interplanar distance in graphite |44j . 



This model Hamiltonian describes with good accuracy 
the lattice parameters and compressibility of graphite 
and diamond at equilibrium, but drastically overesti- 
mates the compressibility of graphite at high pressure. 
As a consequence the theoretical transition pressure to 
diamond is as high as 129 GPa and the volume jump at 
the transition pressure is very small, the volume being 
4.70 A'^/atom for graphite and 4.60 A'^/atom for dia- 
mond (cfr. the ab-initio EOS of graphite and diamond 
of Ref. 0). Nevertheless this model graphite represents 
a good testing case for the new simulation scheme. In 
fact, the small volume jump at the theoretical transition 
pressure prevents the reduction of the activation barrier 
by overpressurization and consequently the transition to 
diamond does not take place in a Parrinello-Rahman sim- 
ulation (70 ps long) even by increasing the pressure up 
to 700 GPa and temperature up to 1000 K. Conversely, 
within the new simulation scheme the transformation to 
diamond occurs very close to the theoretical transition 
density. We have performed a metadynamics simulation 
of graphite at 15 GPa and 300 K with a supercell contain- 
ing 128 atoms initially arranged in the graphite structure 
with four graphitic planes per cell in the ABAB (hexago- 
nal) stacking During the simulation run, 28 ps long, 
we have observed several (of the order of 15) forward and 
backward transitions between graphite and diamond. A 
clear monitoring of the phase transition is given by the 
time evolution of the indicator x which provides a sharp 
distinction between the hybridization states of carbon 
atoms as 



X 



-E-E 



where 



j>k 



E^ 



(10) 



(11) 
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with d = 1.8 A and A = 0.05 A. The index I runs over 
all the N atoms of the simulation cell and is the co- 
ordination number of atom i-th. The index k and j run 
over atoms neighboring to atom i-th and 9jik is the an- 
gle subtended by the tern jik whose contribution to x 
is weighted by the product of the partial coordination 
number riik and n.y. The cosine function in (|10|l is able 
to discriminate between the hybridization sp^ and sp'^. 
In fact the indicator x has a value of -0.125 for graphite 
and -0.06 for diamond at 15 GPa and 300 K. For the sake 
of clarity the time evolution of x is reported in Fig. [7| 
only for the first part of the run. Two transitions from 
graphite to diamond are clearly identified at 10.5 ps and 
at 14.0 ps. Similarly, other twelve transformations can be 
identified in the other 12 ps of simulation. The Gaussian 
potential fills first the free-energy basin corresponding to 
graphite and, once the system is driven to the new phase, 
the basin corresponding to the final structure is then pro- 
gressively filled. Once both basins are filled, the system 
is able to oscillate from one structure to the other. In 
agreement with the ab-initio ParrincUo-Rahman simula- 
tion of Ref.^, some of the structures observed during 
the oscillations from graphite to diamond are a mixture 
of cubic and hexagonal diamond. The transformation at 
10.5 ps (14.0 ps) starts at a volume of 4.45 A^/atom (4.70 
A'^/atom), very close to the theoretical transition volume 
of 4.70 A'^/atom obtained from the theoretical equation 
of state of graphite and diamond. These results demon- 
strate the effectiveness of the metadynamics scheme pre- 
sented above. The phase transition occurs spontaneously 
at the theoretical transition density whereas it does not 
take place in a Parrinello- Rahman simulation (within the 
tight-binding model) even by overpressurizing the system 
up to five times the theoretical transition pressure. This 
metadynamics scheme would be particularly suitable and 
probably superior to the one based on the h matrix (as 
presented in section|nl) for the simulation of phase transi- 
tions described in terms of solid state chemical reactions. 



IV. CONCLUSIONS AND OUTLOOK 

The methodology based on exploration of the free en- 
ergy surface results in a substantial improvement of simu- 
lation of structural transitions in crystals. Using several 
inorganic and organic crystals as examples we demon- 
strated its general applicability as well as heuristic value 
in various cases where the plain constant-pressure sim- 
ulations fail. The main advantages that emerge can be 



summarized as follows. Starting from an initial structure, 
the new approach is able to find a number of realistic 
polymorphs, as demonstrated by the examples of zeolite 
and benzene crystals. For the search for new structures, 
a classical force-field can be used resulting in a simulation 
cost comparable to that of a standard classical MD sim- 
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FIG. 7: Graphite: time evolution of the indicator x (see 
Eq. HIO^ in the text) which discriminates between the struc- 
tures of graphite and diamond. Only the first part of the 
simulation is reported for sake of clarity. Along the whole 
run 28 ps long, we have seen 15 jumps of x similar to those 
reported in the figure which equally correspond to oscillations 
between graphite and diamond. 

phases can be understood at conditions close to those of 
experiment, as shown by the example of graphite. 

Clearly, there still remains a lot of room for improve- 
ment, including the algorithm itself. As we see it, the 
major difficulty in the discrete version of the algorithm 
described in section ^ is the choice of parameters W and 
6h. One of possible alternatives is to use also in this case 
a continuous version of metadynamics with an adaptive 
mechanism for the choice of Gaussian parameters (includ- 
ing anisotropic Gaussians), based on monitoring the cell 
fiuctuations and estimating the shape of the well. The 
work in this direction is in progress. Also, an accurate 
calculation of a transition pressure at finite temperature 
still remains a challenge. The ultimate goal in this field 
should be a crystal structure prediction from "scratch" , 
based only on the knowledge of the chemical composition 
of the system. We believe that methods described in this 
paper represent a substantial step forward towards this 
goal. 
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